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Abstract. In this paper we demonstrate that the numerical method of 
steepest descent fails when applied in a straight forward fashion to the 
most commonly occurring highly oscillatory integrals in scattering theory. 
Through a polar change of variables, however, the integral can be brought 
on a form that can be solved efficiently using a mix of oscillatory integration 
techniques and classical quadrature. The approach is described in detail and 
demonstrated numerically on integration problems taken from applications. 



1. Introduction 

Efficient numerical approximation of oscillatory integrals is a challenging 
problem with a wide range of applications. Multivariate oscillatory integra- 
tion is particularly difficult, with methods and theoretical results in practice 
being limited in applicability to certain classes of integrals. This work concerns 
such a class, namely integrals of the form 

/(x)e''^f{x)dx, 

' D 

where (7(x) behaves locally around a point xq G .D C M" like |x — xo|", and uj 
may be large compared to D. This particular kind of integral appears frequently 
in literature, typically with n = 2,3 and a = 1, e.g. in acoustics [231 E], and 
scattering theory |2H [25l |7] , and more specifically in high frequency BEM [6] . 
The Greens function of the Helmholtz equation possesses this kind of special 
point, and this is clearly one reason for the ubiquitousness of these kinds of 
integrals; any computation of a wave field in space from a surface field, like in 
e.g. BEM, involves integrating the Green's function. 

An extensive amount of research has been devoted to methods for efficient 
evaluation of univariate integrals, integrals of the form, 

rb 

(1.1) / fixy'^six^x, 



where cj is a parameter that may be large. Recently developed methods for such 
integrals include variants of the Filon-type methods [l6l|T5l[ini[l], the numerical 
method of steepest descent [131 12] , and Levin- type methods [201 . For an 



oscillatory integral of the form (1.1), certain special points, namely endpoints, 



singularities, and stationary points, i.e., points where g'{x) = 0, determine the 
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asymptotic behaviour of the integral as a; — t- oo. The above mentioned methods 
treat these asymptotically contributing points specifically such that terms in the 
asymptotic expansion of the error vanish, leading to methods with an error of 
size O{io~^), for some /? > 0. For example, a Filon-type approximation is 
obtained by interpolating f{x) with a number of derivatives at all contributing 
points, and the approximation is constructed a liner combination of moments. 
In this case the asymptotic order /3 of the approximation depends on how many 
derivatives that are interpolated |16j . 

Multivariate oscillatory integrals have the general form 



where now D C M"", and f,g : — )■ M. In this case we can think of the end- 
points in the univariate case as corresponding to corners of dD in n dimensions. 
The notion of a stationary points naturally extends to higher dimensions, being 
a point where 'Vg{x) = 0. In addition, a new sort of contributing point without 
a straight-forward parallel in the univariate case will appear, namely the res- 
onance point, which is a point where 'Vg{x.) is perpendicular to the boundary. 
This, and other added complications, are reason why general-purpose methods 
for multivariate oscillatory integrals can be considered a hard problem. The 
Filon-type methods do generalise in a quite natural way to higher dimensions 
|18l I17j . however, moments are not generally available in higher dimensions, 
and supplementary methods are needed to obtain these. Levin-type methods 
do not rely on moments, but no general procedure for handling cases with sta- 
tionary points and resonance points is available for these methods [23]. The 
numerical method of steepest descent can be generalised to multivariate cases 
by regarding multivariate integrals as nested univariate integrals |14j . and this 
approach can handle stationary points and resonance points to a certain degree. 
This is likely the most general method among the above mentioned. 

In this paper it will however be demonstrated that whenever ^(x) behaves 
locally around a point xq E -D like |x— xo|, this induces problems for the method 
of steepest descent. As such, the only method which seem applicable is the so- 
called localised method of stationary phase due to Ganesh et al. [11], which was 
developed specifically for use in a BEM solver. Here it is suggested that a polar 
change of variables is needed to cancel a singularity of the amplitude function 
/(x) at xq, as is commonly seen in scattering integrals, and the remaining part 
is integrated with a combination of oscillatory and non-oscillatory techniques. 
In this work it will however become apparent that this change of variables is 
also precisely what is needed in order to apply the numerical method of steepest 
descent. 

This paper is built up as follows: Some relevant facts about the numerical 
steepest descent method is presented in section ^ and we shall see demon- 
strated by examples the failure of numerical steepest descent when applied to 
the integrals of our interest, as well as an outline of how to modify the method 
such that it works. In section ^ we shall put these observations on a firmer 
foundation and derive some results that will aid in constructing efficient meth- 
ods for these integrals based in the method of steepest descent. Finally, in 
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Q we demonstrate the resulting methods on one artificial example, and two 
examples taken directly from applications. 

2. Numerical steepest descent and cubature 



Applying the method of steepest descent to the univariate integral (1.1), 
assuming that /(x) and g{x) can be analytically extended to the complex plane, 
and that g{x) is monotone, yields |13] 



/(x)e''^^'(^)dx = G(a)-G(6), 

J a 
where 

roc 

(2.1) G{x) = / f{h,ip))e-'^Ph',{p)dp, 



JO 

where in turn hx{p) is the path of steepest descent, implicitly defined through 
the equation. 

9{hx{p)) = g{x)+ip. 
Note that the requirement that f{x) and g{x) are analytic can be relaxed to 
analytic in appropriate regions of the complex plane [12] . More generally a 
singularity in the complex plane might introduce an exponentially subdominant 
term [13j, and 

/ fixy^'si^Ux = G{a) - Gib) + ©(e-^). 

J a 

In the following the notation 0{e~'^) will be used to indicate any exponentialy 
subdominant term that might be present, depending on the analytic properties 
of the integrand. 

For small p, the path of steepest descent will behave like 



X 



i.e., some brach of this expression, where a — 1 is the number of vanishing 
derivatives of g{x) at the point x [2]. An important consequence of this is that 
a p~''"~^^/"-type singularity is introduced in the steepest descent integral ( |2.1[ ). 
A substitution p — )• will produce an analytic integrand, given that f{x) and 
g{x) are analytic, 



oo 



(2.2) G{x) = a / /(/i.(g"))e"-'? Kiq^q'^'dq 

Jo 

foo fa fa 

= - / f{KC-)K{-)e~h-'''dt. 

W Jq U! CO 

The last integral is on a form that can numerically be resolved efficiently by 
Gaussian quadrature. Given a quadrature rule with m points and weights 
{x",w'j}'j^i which is Gaussian with respect to the weight e~^ on [0,oo), 
i.e., Gauss-Legendre for a = 1, and half-space Gauss-Hermite for a = 2. By 



Lemma 1 from [8j the error of this approximation applied to the integral (2.2) 
is ©(cj-^^'^-i)/"). 
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2.1. Cubature. The method of steepest descent is inherently a procedure for 
univariate integrals; no higher dimensional extension of the paths of steepest 
descent are known. The method can however be developed into a cubature rule 
by regarding a double integral as nested univariate integration. Huybrechs & 
Vandewalle in [T3] demonstrate this approach on cases with stationary points, 
resonance points and corner points, giving a method with high asymptotic order. 

In the following example we shall however see that the numerical method 
fails on the particular type of integral that we are interested in. 

Example 1. From [14J we have 

[ [ f{x,y)e'^'^^^dydx = F{0,0)-F{a,0)-F{0,b) + F{a,b) + O{e-'"^). 
Jo Jo 

where, 

(2.3) 

F{x,y)= / /(n,,,(p),z;,K,,(p),g))<,^(p)-^(n,(p),g)e-(^'+'?)dpdg. 
JO Jo (^1 

The paths of steepest descent for the inner integral are easily obtained, and we 

get 

Vy{x,q) = \Jy'^ - + 2\qyJ x'^ + 

Ux,y{v) = \/ x"^ - p"^ + 2ipyGp~^y^. 

Now, assume f{x, y) = 1 and let us then study the contribution from the origin 
in particular, 

F(0,0) ^_ r f"^ P + Q e-'^iP+g)dpdq. 
Jo Jo V 2pq + g2 

It can be showed that F{0, 0) = There is a sqare-root singularity at q = 0, 

which follows from the fact that the oscillator function has a line of 

stationary points along y = 0. Substituting q ^ (f' , will resolve this, 

poo roo »5 I «2 

F(0,0) = -2/ / ^ ^ ^ e-^(P+^ )dpdg 
Jo Jo V 2p + g2 

roo roo + \ +2 

= -2u-^ / 4i±i^e-(*^+*i)dtidt2. 
Jo Jo yj2ti + tl 

Now it's apparent that this method can never have asymptotic order higher than 
2 unless this integral is resolved exactly by the quadrature method. In other 
words, the relative error will generally not decrease with uj. Moreover, a pole 
is present at q = ±.\y/2p, which will also reduce the efficiency of any standard 
quadrature quadrature rule, like e.g. Gauss-quadrature. 

The failure of the steepest descent method in this case should not come as a 
surprise, since the oscillator function ^/x"^ + y'^ is not differentiable at the origin. 
Indeed, the line of stationary points on the y-axis actually degenerates into a 
regular endpoint as x — )• 0. Now, the observant reader has probably already 
realised that a polar change of variables would swiftly resolve the problem in 
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this example. One way of seeing this is by observing that, in the sense of a 
regularised integral, 

F(0,0)=/ / ei^v^'+2''dydx= / / e^'^Vdrd^, 
Jo Jo Jo Jo 

which is, of course, the way one arrives at the exact value of However, 

this form also turns out to be much better suited for computations. Let us 

remove the assumption that f{x,y) = 1. Now, again in sense of a regularised 

integral, 

/•oo roo . /"^/^ /"OO 

F(0,0)=/ / /(x,y)e^^v^'+2''dxdy = / / f{r,9)e^'^Wdrde, 
Jo Jo Jo Jo 

where 

f{r,6) := f{r cos 6, r sinO). 

The method of steepest descent applied to the inner integral gives, assuming 
that /(r, 9) is analytic in a sufficiently large part of the complex plane in its 
first argument, 

/'7r/2 rco 

F(0, 0) = i / e^'^Pfiip, 9)dpd9. 
Jo Jo 

Note that the outer integral is non-oscillatory as a function of 9. Thus, the outer 
integral can be resolved with a classical quadrature rule. The inner integral is 
efficiently handled by scaled Gauss-Laguerre quadrature. In other words, the 
polar change of variables has effectively confined the oscillatory behaviour to 
the radial dimension. In the following we shall see that this is not limited to 
this particular oscillator function. Also note that if f(x,y) ^ (x^ + y'^)~'^^'^ 
near the origin, then f{r,9) ~ r~'^'^^, thus cancelling the 1/r-type singularity 
typically seen in scattering integrals. Finally note that the integration range 
of the ^-integration reflects the shape of the domain around the origin. If the 
origin was contained in the domain, then the outer integration range would be 
[0,27r], and we would be integrating a periodic function, which means that a 
simple trapezoidal rule would be the natural choice of quadrature in this case. 

3. Theoretical results 

The method outlined in the example in the previous section works under the 
premise that the contribution from the special point xo can be separated, and 
that it is of the form 

rf} roo 

/ / f{r,9)e''^3^'''^^drd9. 
Jo Jo 

In this section we shall prove some fairly general results regarding the appli- 
cation of the method of steepest descent in conjunction with a polar change 
of variables. The results are formulated as a kind of pre- quadrature rule, by 
which we mean that only the problematic part has been treated. No choice 
of quadrature rule has been made for the rest, but existing methods are here 
applicable. Though most real-world cases are with n = 2,3, the results are 
formulated for integration in M", since this implies little extra work. The tool 
here is n-spherical coordinates |[3], which is the higher dimensional equivalent 
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of the well-known polar and spherical coordinates. Here the substitution takes 
the form, 

i-i 

Xj rcos{ipj)Y[sm{ipi), j = l,...,n- 1, 
1=1 

n—l 

rY[sin{(pi), 
1=1 

where (pj G [0, tt], j = 1, ... ,n — 2, and tpn-i € [0,27r]. In vectorial form this 
can be written 

where S"~^ denotes the (n — l)-sphere. The volume element has the simple 
form r"~^dG, where, 

n-l 

1=1 

In the following, we shall assume that the special point xq = 0, and that 

5f(xo) = 0. Note that this implies no lack of generality, since this represents 
a simple translation of the coordinate system, plus a scaling of the integral: 
/ /(x) exp(iwc/(x))da; = exp(ia;5r(xo)) / /(xq + x) exp(iw(c/(xo + x.) - g{^o))dx. 

3.1. Unbounded integrals. For the unbounded case we consider a domain of 
integration that is an infinite cone defined in terms of a subset of the (n — 1)- 
sphere, W C S"~^. Note that this case includes the case of the entire Euclidean 
space as a special case. 

Lemma 1. For n > 1, assume that the integral, 

nf] ■■= f /(x)e-^(-)dx, 
Jd 

exists, with D being an infinite cone: D = {rO, where O G C S""^, r G 
[0, oo]}. Suppose that for the oscillator function ^(x) the following conditions 
hold whenever x G D/{0}: 

(1) g{x) is continuously differentiable. 

(2) V<7(x) / 0. 

(3) l^g{r@) > 0. 

Assume in addition that r"'~^f{r@) and g{r@), @ G W, are analytic functions 
of r in a neighbourhood of the origin. Then 

(3.1) /[/] = - / r fipoip, e)@)^ip, e)e-^Pdpde + o(e— ), 

n- Jw Jo 

where po{jp,0) satisfies the equation 

Proof. Writing the integral in n-sphcrical coordinates, we get, 



/[/]= / /(re)e^'^^('^®V"-Mrde. 

Jw Jo 
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Let US apply the method of steepest descent to the inner integral. Now condi- 
tions 1, 2, and 3 on (?(x) ensures that there are no contributing points of g{rQ) 
for r > 0, and condition 3 ensures that po{p, Q) exists for all 0, thus, 

-R 1 roo r)( rP\ 

n Jq dp 



n 



fiPRip, 9)9) ^(p, 9)e-'^^dp + 0{e-n 



where g{pR{p,&)&) = g{R@) + ip- Letting i? — )• oo the second integral must 
vanish. The result follows. □ 

From this Lemma a very simple, and potentially highly efficient, quadra- 
ture rule, or more precisely a pre-quadrature rule, for these particular kinds of 
integrals follows. 

Theorem 1. Assume the conditions of Lemma^ are satisfied. Let a be the 
smallest integer such that 

d^g{rQ) 
and assume that 

C*r" r=0+ 

Let {x",w'j}J^i be the weights and nodes of the m-point rule being Gaussian 
with respect to the weight e~^" on [0,oo). Define, 

(3.4) Q,[/](9) = —f2w,xJ-'fipoix'^/u;,e)e)^{xJ/u;,Q). 



(3.2) 



= 0. V9gVFcS"~\ / = 0,l,...,a-l, 

T- = 0+ 



. , a"c/(r9) 
(3.3) ^ 



Then 



nu! ^-^ JO Qp 

m - [ Qr[fmde = o{u'"-^). 

Jw 



Proof. Departing from Eq. (3.1), we consider the inner integral 

hnner{@) := - H f{po{p,e)e)^{p,e)e~^Pdp. 

n Jq dp 

From of conditions (3.2) and (3.3) it follows that p{p,Q) has an expansion of 
the form 

po{p, 9) ~ ai(9)pi/" + 02(9)//- + . . . , 
valid for all 9 G VF. Introducing the change of variables p ^ q°' we get 



linneriQ) = " / /(po(g", 9)9)^((?°, 9)(z"-^e-"^ dq, 

n Jq dp 
whose integrand is analytic around q = 0. Now, applying Lemma 1 of [S], 

I^nneriQ) - — W^xf'^ f (poix'^ / , 9)9)^ (x'^ / U , O) = ©(o;"^). 

Integrating over W yields the sought conclusion. □ 
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Remark 1. Note here that ~ q"'~ , so in the case where /(x) is regular at 

the origin, applying a Gaussian rule with the weight x^~^e~^" will yield higher 
asymptotic accuracy. 

The following corollary follows directly from the above theorem, and im- 
plies that when a quadratm'e method is applied for the outer integration, the 
error will decrease asymptotically with u down to the error level of the (non- 
oscillatory) outer integration. 

Corollary 1. Let Qw ■ W be a quadrature rule, 

h{e)de = Qw[h] + E[h]. 

w 
Then 

nn - Qw[Qr[f]] = E[Qr[f]] + o{u-"-^). 

3.2. Bounded integrals. Next we shall investigate the case of integration over 
a more general bounded domain D C M", 



/[/] :=^/(x)e-3Wdx. 



'D 

Assume that D is star shaped with respect to the origin, i.e., the segment 
(0, x) is contained in D for all x G dD, and we can write D = {rQ, Vr G 

In n-spherical coordinates centred at the origin we then have, 

I[f]= // / /(re)e''^^'('^®V"-Mrde. 

JJw Jo 

Assuming the functions r"~^/(r0) and g{r@) are analytic as a function of r in 
a neighbourhood around the origin, just as in the proof of Lemma [T| we can 
apply the method of steepest descent and get 

(3.5) m = -f[ r f{po{p,@)e)^ip,e)e-'^PdpdQ 

n JJw Jo op 

-- [[ e'-^(®) r f(p^(p^e)e)^i^(p,e)e-Pdpde + o{e-n, 

n JJw Jo op 

where pnip, 0) is defined through the equation g{pR{p, ©)) = R{@) + ip. This 
leads to a theorem of the same kind as Theorem [T| but with an additional term 
that includes the contributions from the boundary. 

Theorem 2. Assume D is of the form D = {tQ, Vt G [0,R{Q)], Q £ W C 
S"-^}, Ve G ly C S"-^ Furthermore assume that the conditions 1, 2, and 3 
in Lemma^ hold for X. G D/{0}, and that conditions (3.2) and (3.3) hold. Let 
{x'j , w'j}'JLi be the weights and nodes of the m-point rule being Gaussian with 
respect to the weight e~^" on [0,oo) and let Qr[f] be given by equation (3.4). 
Now define 

(3.6) Qrifm = Qrlfm - ^e-«(®)^"/(p«(p,G)e)^(p,e)e-fdp. 
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then 

-{2m~l)/a\ 



[[ QrUmde - I[f] = 0{u;' 
JJW 



Proof. Under the given assumptions, we have the that Eq. (3.5) holds, and we 



simply apply Theorem [T| and the result follows. □ 

Let us end this section with some remarks 

Remark 2. Only for the central contribution will we see the oscillations being 
confined to the radial direction, and in this case must the remaining integral 
be treated with oscillatory quadrature techniques with sufficiently high order in 
order to retain the asymptotic error decay indicated in here. In practice, for 
the method of steepest descent to be directly applicable to the remainder we need 
that R{Q) > 0, V0 G W . That means cases where the origin is contained in 
the interior of D, and intersections of such sets with cones, as are discussed in 
Theorem^ Note that this stronger condition does not allow the origin to be on 
a curved section of the boundary. 

Remark 3. Theorem^ gives a notably different decomposition that what was 
seen in Example [7| In the example the method of steepest descent is applied to 
the integral given in it's original coordinates, and then the central contribution 
is replaced by an unbounded integral on polar form. However, if the method 



of steepest descent is applied to the remainder in (3.6), this will result in a 
decomposition like that of Example [I[ with contributions being local to each 
corner of the rectangle. The contributions will be given in different coordinates, 
but at least asymptotically the decompositions must be the same. This will be 
touched upon in the following section. 

4. Numerical experiments 

Let us first consider an artificial example. The following integral has a rela- 
tively simple closed form solution, 

/ / ^dzdydx, 

-oo i-oo i-oo (x2 + 2y2 + 3z2)(i + y^2^2FT3^) 

-7r(i cos(a;) + sin(tj))(7r + 2iCi(a;) — 2Si(a;)). 

The spherical substitution is x — )• r cos c/^i, y — t- r sin ipi cos ip2, z ^ r sin ipi sin ip2, 
never mind that a scaling of each dimension would in this case produce a simpler 
expression. Now, in spherical variables. 



g{r, (pi,(f2) = r\l cos'^{(pi) + ^(5 + cos(2932)) sin^{(pi] 



2 

ip 



and then 

p{p,ipi,(p2) = 

Jcos'^{(fi) + ^(5 + cos(2v32)) sm'^{(fi) 

Theorem [T] indicates that a method with high asymptotic order is obtained 
by applying a Gaussian rule in the radial direction. For the given integral we 
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Figure 1. Absolute errors in resolving the integral (4.1). 
Straight lines indicate lines of ~ w^^, ~ w^^, ~ uj^^'^ and 
~ u)~^'^ , respectively. Dashed lines indicate result of same ex- 
periment with = 30 points in each dimension for the outer 
integration. 



apply formula (3.4), where we use Clenshaw-Curtis quadrature with 50 points 
in the c/^i-direction, and the trapezoidal rule with 50 points in the (periodic) 
(/?2-direction for integrating over the sphere. The function that here is being 
integrated is visualised in Fig. [2j The result of this experiment can be seen in 
Fig. [T| where the method has been applied with the number of complex Gauss- 
ian points in the radial integration being between 2 and 8. One observes an 
asymptotic decrease in the error matches the prediction of Theorem [TJ Observe 
that machine precision can be attained with relatively few quadrature points 
when u! is sufficiently large. In Fig. [T] curves are also included that show the 
result of using less points (30) for the outer integration. These curves match 
what should be expected from Corollary [TJ 



4.1. A problem from acoustics. For the problem of a rectangular duct ter- 
minating in an infinite baffle, sound pressure at the end of the duct expressed 
in terms of the axial velocity of the baffle, v{x,y) is of the form (see [19] for 
details), 

p{xo,yo) = C / , v{x,y)dydx, 

Jo Jo y/ [X - xo)^ + [y - yoy 
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-s 0.5^ 



CO 




Figure 2. The absolute value of ImneriO-, ^)-, for k = 100. 




Figure 3. Splitting of the outer integration. 



for some constant C . We shall consider a model case of this integral, 

fa ^-b gitj-y/ x^+y'^ 



JO 



^/x"^ + y' 



;f{x,y)dydx, 



For simplicity we denote /(r, 0) = f (r cos{9) , r sm(9)) . The central contribution 
is now computed from (3.4) 



Qr[f] = -^wJC^,9), 
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and an approximation to the full integral is, by By Theorem [2j 

Qr[f]de = r'^ Qr[f]de - hAfl 

Jo Jo 

where 

/'7r/2 roo 

leAf] = i / e'^^^W / f{R{0) + ip, 9)e-^Pdpde. 
Jo Jo 

In order to compute the last integral, we apply the method of steepest descent 
again, noting that the integral should be split in the 6 direction. Let in the 
following /3 = arctan^, as illustrated in Figjsj then Iext[f] = I^xtif] + ^extif] 
with 

Ilxtif] = i / e-'^-^(^) / f{asec{d) + ip,e)e-'^Pdpde, 
Jo Jo 

and 

/'7r/2 roo 

Iext[f] = '^ gi<^f.csc(e) / f(^t)csc{e)+[p,e)e-^Pdpde. 

J 13 Jo 

We compute two paths of steepest descent for the ^-integration in each of the 
integrals, 

^1,1(9) = sec"^(l + iq/a), hio(q) = sec"^ ( ^ 

V a 

— - — j, /ii,2(g) = csc"^(l + ig/6). 
where rj := a sec (3 = h esc (3 = \/ + 6^ . Thus 

louter [f] = I outer [/] + outer [/] ~ outer [/] ~ ^ outer [/] + outer [/] ~ ^ outer [/] ' 

where 

(a + i(7) Y^2iga — ^2 



/•oo .oo/(r? + ig + ipsec-M2±i2 ) 

llx'tif] = -ae'-' / ' L ^e-(^+^)dpdg, 

Jo [i] + iq)y'b'^ — q"^ + 2iqr] 

42^1 [/] = te'"'' / / , / ^ -e-"(P+g)dpdg, 

JO (r/ + ig) y - g2 _|_ 2iqr/ 

Jo Jo (ft + ig)V2ig6 - 

Now one should note that the integrals lifter if] ^^^^ ^outer if] have a square-root 
singularity at g = 0, this is because the points [a, 0] and [0, b] are resonance 
points, which induce stationary points in R{9). Substituting p ^ p/uj and 
q — )• (f' /oj in these integrals, and p ^ pjuj and q — )• qjuj in the other two 
integrals bring them on a form that can efficiently be evaluated with tensor- 
product half-space Gauss-Hermite x Gauss Laguerre, and Gauss-Laguerre x 
Gauss Laguerre. 
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To test the method we apply it to the following problem, 

/ / y cos{x)dydx = - / (e'^^ - e''^^'^+^) cos(x)dx. 

Jo Jo y/x^ + y"^ ^ Jo 

By integrating the right hand integral to machine precision with Matlah 's quadgk 
routine we have a reliable reference solution. 

In Figure |4] we observe that on this problem machine precision can be reached 
with relatively few points in the complex plane. For producing this figure 
the non-oscillatory integration was performed with 30-point Clenshaw-Curtis 
quadrature. Note that using the decomposition from Example [T] with modified 
central contribution (dashed lines) gives a very similar picture. 




5 10 20 50 100 



CO 

Figure 4. Relative error using ugl = 2, 4, 6, 8 Gauss-Laguerre 
points and hgh = 2nGL half-space Gauss-Hermite points. Dot- 
ted lines indicate a direct application of numerical steepest de- 
scent (cf. Ex. [T]), and dashed lines idem, but with modified 
central contribution. 

4.2. A problem from high frequency scattering. The scalar exterior scat- 
tering problem is often formulated in terms of the Kirchhoff-Helmholtz integral 
equation for an unknown surface field (see e.g. [Hj for details). Here we shall 
limit the discussion to the case of Dirichlet boundary conditions, where one can 
use a Fredholm equation of the first kind, 

/ g(x)G(x-y)dSx = -n^(y), y £ dn, 
Jan 

where d^l denotes the surface of the scatterer, G{x, y) the free-space Green's 
function, ti*(y) is the incident field, and q{x) is the unknown surface field. We 
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shall in the following discuss the 3D case with an incident plane wave, 

pifc|x-y| 

G(x,y) = — 

47r|x — y| 

u\y) = e'''^y'^\ where ||d|| = 1. 

The field q{x.) is in general oscillatory with wavelength k, however, as has 
been demonstrated in e.g. [9], by factoring out the phase of the Kirchhoff 
approximation, q{'y) = q{:x.)e^''^''^ , the remaining amplitude function q'(x) will 
be non-oscillatory in illuminated regions of the scatterer. The new unknown 
satisfies 

giA:{|x-y| + (x-y)-d) 

g(x) — dsx = -1, y £ dn. 

lan 47r|x-y| 

Now, relatively few degrees of freedom are needed to represent Q'(x), but a 
discretisation matrix will necessarily involve integrals with the same kernel as 
the integral here on the left. 

Studying the case of a spherical scatterer of radius 1 centred at the origin, 
we set X = [1, 0, 0]^, and investigate the resulting integral for different incident 
directions d = [di, d2, d^] . The corresponding integral equation is 

(4.2) 

r>7r /•2tt ife(-y/2+2cos(<)!)2) sin(?ii)— sin(<^i) cos(</)2)cii— sin((/)i) sin(<7i2)rf2 — (l+cos(0i))d3) 
■H>2) 



-'0 



4TTy/2 + 2 COs((/)2) sin((/)i) 

sm{<pi)d<f>2d-(f>i = -1. 



Using the polar coordinate substitution, centred around 01 = 7r/2 02 = 0, on 
this integral we can show that 



dg{rQ) 



i + [-d3,d2fe, 



lr=0+ 

and it follows from this that the conditions of Theorem [2] are satisfied whenever 
di 7^ 0, which translates to x being on the shadow boundary. Computing only 
the central contribution, i.e., Qr[4>\dO, gives a local approximation of the in- 
tegral, much like what is achieved with the localised method of stationary phase 
. A way to test the usefulness of this approximation is to consider a 3D ver- 
sion of the method in [3], whereby a prototype asymptotic Filon-approximation 
of q would be given as, 

g(x) . 

Here wq is constructed as a local approximation of the integral of the left hand 
side of (|42]) with g(0i,02) = 1- In Table [l] we see the result of such a test, 



where now d = [cos(^'), 0, sin(^')]^, and we compute the relative error w.r.t 
the analytic solution given by the Mie series |21j . For this table the number of 
Laguerre points are m = 5, and the angular direction is computed with a 100 
point trapezoidal rule. This test clearly shows the potential of steepest descent 
for high feequency scattering computations. 



STEEPEST DESCENT APPLIED TO SCATTERING INTEGRALS 

Table 1. Relative error of prototype Filon-approximation to q'(x). 
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10 

TT 

5 

TT 

3 



50 

5.43380e-04 
2.28316e-03 
1.227586-02 
4.79852e-02 



100 
1.37074e-04 
1.11796e-03 

6.40899e-03 
4.05070e-02 



150 
6.10265e-05 
7.42382e-04 

4.31359e-03 
3.24495e-02 



200 
3.43482e-05 
5.56018e-04 
3.246776-03 
2.678786-02 
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